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Abstract 

Reaumuria soongorica (Pall.) Maxim., a resurrection semi-shrub, is a typical constructive and dominant species in desert 
ecosystems in northwestern China. However, the gene expression characteristics of R. soongorica under drought stress have 
not been elucidated. Digital gene expression analysis was performed using lllumina technique to investigate differentially 
expressed genes (DEGs) between control and PEG-treated samples of R. soongorica. A total of 212,338 and 211,052 distinct 
tags were detected in the control and PEG-treated libraries, respectively. A total of 1,325 genes were identified as DEGs, 379 
(28.6%) of which were up-regulated and 946 (71.4%) were down-regulated in response to drought stress. Functional 
annotation analysis identified numerous drought-inducible genes with various functions in response to drought stress. A 
number of regulatory proteins, functional proteins, and proteins induced by other stress factors in R. soongorica were 
identified. Alteration in the regulatory proteins (transcription factors and protein kinase) may be involved in signal 
transduction. Functional proteins, including flavonoid biosynthetic proteins, late embryogenesis abundant (LEA) proteins, 
small heat shock proteins (sHSP), and aquaporin and proline transporter may play protective roles in response to drought 
stress. Flavonoids, LEA proteins and sHSP function as reactive oxygen species scavenger or molecular chaperone. Aquaporin 
and proline transporters regulate the distribution of water and proline throughout the whole plant. The tolerance ability of 
R. soongorica may be gained through effective signal transduction and enhanced protection of functional proteins to 
reestablish cellular homeostasis. DEGs obtained in this study may provide useful insights to help further understand the 
drought-tolerant mechanism of R. soongorica. 
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Introduction 

Water deficit is one of the most significant abiotic stresses that 
influence germination, growth, development, and productivity of 
plants [1]. Drought stress causes stomatal closure, limited gas 
exchange, and reduced photosynthesis in plants [2]. Over- 
reduction of photosynthetic electron transport chain induces tlie 
generation of reactive oxygen species (ROS), such as singlet 
oxygen ('O2), superoxide anion (O2 ), hydrogen peroxide (H2O2), 
and hydroxyl radical ('OH), which damage cellular structures and 
macromolecules [3]. Plants have enzymatic and non-enzymatic 
systems to eliminate ROS. Moreover, plants have developed 
drought-resistance strategies such as succulent leaves, formation of 
osmophUic globules, stomatal movement, and reduction of leaf 
water potential [4,5]. The understanding of plant responses to 
water deficit have improved because of the application of 
molecular techniques. The expression of numerous genes in 
response to drought stress was described in previous studies [6,7]. 
Seki et al. [8] classified the genes expressed during stress into two 



groups: (i) genes encoding proteins involved in signal transduction 
(protein kinases and transcription factors) and (ii) genes with 
products, such as late embryogenesis abundant (LEA) proteins, 
chaperone, osmoprotectants, and detoxification enzymes, that 
directly protect cells against stress. Genomic and transcriptomic 
analyses revealed that various transcriptional regulatory systems 
are involved in stress-responsive gene induction. Several different 
sets of cis- and ten,r-acting factors are known to be induced by 
drought stress at the molecular level [9]. A large number of 
metabolites and proteins have also been reported to be up 
regulated in response to drought stress [10]. 

Reaumuria soongorica (Pall.) Maxim, is an extreme xerophytic 
semi-shrub and a typical constructive and dominant species in the 
desert vegetation in China. R. soongorica forms the zonal landscape 
and is widely distributed in northwest China [1 1]. Water supply is 
one of the main limiting factors in the habitat of this species. The 
unique adaptive strategies in the morphology and physiology oi R. 
soongorica, such as thick cuticle, hoUow stomata, and accumulation 
of some low-molecular-weight metabolites, are attributed to its 
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special distribution area [4]. R. soongorica is a vascular flowering 
plant with desiccation-tolerance. R. soongorica leaves wither and 
enter a state of dormancy during dehydration but is revived when 

water becomes available. Thus, R. soongorica is utilized as a valuable 
non-model species in exploring drought-tolerance mechanisms. 
Liu et al. [4] found sucrose, malate, and proline, which play 
important roles in osmoregulation in R. soongorica during water loss. 
Studies on the protection mechanism for photosynthetic proper- 
ties, activity of antioxidant enzyme, and metabolite changes under 
drought stress have also been conducted [12,13]. However, the 
molecular m(;chanism of drought tolerance in R. soongorica is still 
poorly understood. A preliminary genome survey, which included 
genome size, chromosome number, and karyotype, was conducted 
with R. soongorica [14]. In our previous study, the mitogen-activated 
protein kinase gene (RsMPK2) was isolated and its participation as 
possible mediator under different stresses was investigated [15]. 
We also identified and characterized the response of flavonone 3- 
hydrolase gene {RsF3H), a key enzyme-encoding gene involved in 
flavonoid biosynthesis pathway in drought treatment [16]. 
Recently, the transcriptome of R. soongorica was sequenced and 
analyzed to obtain valuable genetic information for the investiga- 
tion on the molecular mechanism of drought tolerance [17]. 
Genetic information from the transcriptome of R. soongorica was 
also used to identify differentially expressed genes (DEGs) under 
drought stress to understand the mechanism of drought tolerance. 

In recent years, with the increasing availability of sequence data, 
expression profiling has been used to identify DEGs in different 
tissues, organs, developmental stages or mutants and determine 
the expression patterns induced by stress in a large number of 
model and non-model organisms. Microarray analysis and 
suppression subtractive hybridization are powerful tools for 
isolating differentially expressed cDNAs that mediate complex 
biological processes in higher eukaryotes [18]. These technologies 
have been used in a number of animals and plants, such as 
Aedesaegypti [19], soybean [20], Populuscanadensis [21], Scylla 
paramamosain [22], cucumber [23], and navel orange [24]. Next- 
generation sequencing technologies offer new approaches for 
global measurements of gene expression with the advantage of 
high efiiciency and low cost. Digital gene expression (DGE) tag 
profihng is based on ultra-high-throughput sequencing of cDNA 
fragments that unir|uely tag the corresponding gene, thereby 
allowing direct quantification of transcript abundance [25] . DGE 
technology allows identification of millions of DEGs without the 
need for prior annotations and permits the analysis of organisms 
that lack genomic information. Thus, DGE has been widely 
utilized to monitor the differences in transcriptional responses 
among different tissues and organs in response to stress in 
silkworm [26], Sagittariatrifolia [27], Dugesia japonica [28], Populus 
[29], cotton [30], spruce [31], Brassica napus [32], and moss [33]. 

In this study, DEGs in leaves oiR. soongorica in response to PEG- 
induced drought stress were examined using DGE tag profiling 
technology. Analysis of gene expression related to stress response 
provides further insight into the molecular mechanisms of stress 
tolerance in R. soongorica. Based on the putative functions of the 
identified genes, some important genes may be cloned. Moreover, 
the cloning of stress tolerance genes and determination of their 
expression patterns under drought stress may reveal attractive 
candidate genes and valuable information to improve drought 
stress tolerance of plants through genetic engineering. 



Materials and Methods 

Ethics Statement 

R. soongorica is widely distributed in northwest China and is not 
listed as endangered or protected species. No specific permissions 
are required for sample collection in the Northern Mountain of 
Lanzhou, and that the field studies did not involve any endangered 
or protected species. 

Plant Materials and Drought Treatment 

S(;cds of R. soongorica were collected from the Northern 
Mountain of Lanzhou, Gansu, China (36°17'N, 103°48'E, 
1700-1900 m elevation). Seeds were planted in the breeding base 
at the foothill. Three-year-old plants were transferred to a 
greenhouse under the following conditions: 150 |J,mol m~^s~' 
(fluorescent tubes), photoperiod 16 h light/8 h darkness, day/ 
night temperature of 25°C/20°C. The plants were washed with 
deionized water to remove surface contaminants. Soil attached to 
roots was also rinsed off with deionized water. The plants were 
then cultivated in 0.125 xMurashige and Skoog medium solution 
(without vitamins and organics) with insufflating air for two weeks 
to accommodate the liquid environment during treatments. The 
samples were classified into untreated and PEG-treated groups. 
The untreated samples were maintained in cultivating conditions 
as previously described. For drought treatment, 15% PEG6000 
was added to the culture solution of R. soongorica. All plants were 
treated for 12 h, from 6:00 to 18:00 local time. Leaves from each 
treatment were collected and immediately stored in liquid nitrogen 
for later use. 

DGE Library Preparation and Sequencing 

Total RNA samples were isolated according to modified CTAB- 
method [34] . Quality and quantity analysis of total RNA, library 
construction, and sequencing were performed at the Beijing 
Genome Institute at Shenzhen, China. An extract of 8 |a,g of total 
RNA was obtained and treated with oligo (dT) magnetic bead 
adsorption to purify mRNA. Oligo (dT) was then used as primer in 
the synthesis of first- and second-strand cDNA. The 5'-ends of tags 
can be generated by two types of endonucleases, namely, Nlalll or 
DpnII. The bead-bound cDNA was subsequently digested with 
restriction enzyme Nlalll, which recognizes and severs CATG 
sites. Fragments other than 3'-cDNA fragments connected to oligo 
(dT) beads were washed away, and the lUumina adaptor 1 was 
ligated to the sticky 5 '-end of the digested bead-bound cDNA 
fragments. The junction of the lUumina adaptor 1 and the CATG 
site constitutes the recognition site of Mmel, which is a type of 
ciidoiiuclease with separate recognition and digestion sites. Mmel 
breaks 1 7 bp downstream of the CATG site, producing tags with 
adaptor 1. After removing 3' fragments with magnetic bead 
precipitation, the lUumina adaptor 2 was ligated to the 3'-ends of 
tags, thereby acquiring tags with different adaptors of both ends to 
form a tag library. After linear PCR amplification, fragments were 
purified by PAGE gel electrophoresis. During the quality control 
steps, Agilent 2100 Bioanalyzer and ABI Step One Plus Real- 
Time PCR System were used in the quantification and qualifica- 
tion of the sample library. Finally, the library was sequenced using 
Illumina HiSec} 2000 sequencer. The available raw data was then 
deposited in the NCBI Gene Expression Omnibus (GEO) 
database (http:/ /www.ncbi.nlm.nih.gov/ geo/query/ acc.cgiPacc 
= GSE55412). 

Analysis and Mapping of DGE Tags 

Sequencing-received raw image data were transformed by base 
calling into sequence data. Prior to mapping the readings from the 
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reference database, all sequences were filtered to remove 3' 
adaptor sequence, low-quality sequences (tags with unknown 
sequences 'N'), tags with a copy number of one (probably 
sequencing error), and tags that were long or too short, leaving 
clean tags of 2 1 nt. For annotation, all clean tags were mapped in 
reference sequences, and a mismatch of only 1 bp was considered. 
We obtained the combined transcriptome of R. soongorica with 
(suffix is "_3") and (sufTrx is "_A") (http://www.ebi.ac.uk/ 
arrayexpress/experiments/E-MTAB-1543/) [17], which was used 
as reference database. Clean tags were mapped in reference 
sequences. The. clean reads mapped in reference database from 
multiple genes were filtered. The remaining clean tags were used 
as unambiguous clean tags. The number of unambiguous clean 
tags for each gene was calculated and normalized to TPM 
(transcripts per million clean tags). 

Identification and Functional Annotation of DEGs 

We used a rigorous algorithm method to identify DEGs 
between two samples. False discovery rate (FDR) was applied to 
determine the threshold of P-value in multiple tests and analysis. 

The DEGs were obtained through FDR <0.001 and |log2Ratio| 
£1. Gene ontology (GO) and pathway classification was used to 
determine the possible functions of all DEGs by mapping the GO 
(http://www.geneontology.org/) and KEGG (http://www. 
genome.jp/kegg/) databases. 

Results 

Analysis of DGE Libraries 

The lUumina platform was used to perform high throughput 
tag-seq analysis on R. soongorica to investigate transcriptome 
responses to PEG-induced drought stress. Total RNA isolated 
from untreated and PEG-treated groups were named as libraries C 
and P, respectively. The major characteristics of the two libraries 
are summarized in Table 1. An average of 7,231,417 total 
sequence tags per library was obtained with 466,957 distinct tag 
sequences from the two libraries. Prior to mapping tag sequences 
to the reference sequences, adaptor tags were filtered (low-quality 
tags and tags with one copy) and produced approximately 
6,969,782 total clean sequence tags per library with 211,695 
distinct clean tag sequences. The distributions of the total and 
distinct clean tag copy numbers showed highly similar tendencies 
in the two libraries (Figure 1). Among the distinct clean tags, more 
than 60% of the distinct clean tags had 2 to 5 copies, 34% of the 
distinct clean tags had 5 to 100 copies, and 4.5% had copy 
numbers higher than 100. 

Analysis of Tag Mapping 

A reference gene database that included 94,878 sequences of the 
R. soongorica unigene was preprocessed for tag mapping. Among 
the sequences, genes with a CATG site accounted for 77.89%. To 
obtain the reference tags, all CATG+17 tags were used as gene 
reference tags. Finally, 218,791 total reference tag sequences with 
182,184 unambiguous reference tags were obtained. Approxi- 
mately 40.44% and 38.75% of the distinct clean tags were mapped 
unambiguously in the unigene database, and 47.94% and 50.2% 
of the distinct clean tags were not mapped to the unigene virtual 
tag database in libraries C and P, respectively (Table 1). 

The sequencing saturation was analyzed in the two libraries to 
estimate whether the sequencing depth was sufficient for the 
transcriptome coverage. The genes that were mapped by all clean 
tags and unambiguous clean tags increased with the total number 
of tags. However, when the sequencing counts reached three 
million tags or higher, the number of detected genes was saturated 



(Figure SI). Given that Dlumina sequencing can distinguish 
transcripts originating from both DNA strands and the strand- 
specific nature of the sequencing tags obtained, we found 
approximately 27% and 26% of distinct tags, which were mapped 
insense transcripts of libraries C and P, respectively (Figure S2). 
Approximately 25% and 24% of distinct tags perfectly matched 
antisense transcripts, which suggested that antisense genes play 
important roles in the transcriptional regulation of drought 
response in R. soongorica. 

Drought-caused Changes in the Global Gene Expression 
in R. Soongorica 

To obtain the transcriptional changes in PEG-treated R. 
soongorica, a rigorous algorithm method was applied to identify 
DEGs from normalized DGE data using pairwise comparison 
between the two groups (C and P). Results showed that 1,325 
genes had FDRs SO. 001, and |log2Ratio| Si, which were 
declared to be the DEGs between the control and PEG-treated 
plants. Among these genes, 379 (28.6%) genes were up-regulated 
and 946 (71.4%) were down-regulated in response to drought 
stress. The expression of most genes was suppressed after drought 
treatment. 

Functional annotation was performed to assign the genes of R. 
soongorica with GO terms. The main GO terms included "cellular 
component", "molecular function", and "biological process". A 
total of 535 DEGs were annotated as "cellular component". The 
categories of "cell part" and "cell" had the same proportion 
(88.0%), followed by "intracellular" and "intracellular part", 
which accounted for 78.5% and 77. 9"^ of assignments, respec- 
tively (Table SI). In addition, a total of 570 DEGs were annotated 
to "molecular function" terms. The categories "catalytic activity", 
"organic cyclic compound binding", and "heterocyclic compound 
binding" were prominent with 60.5%, 25.3% and 25.3% of 
assignments, respectively (Table SI). Furthermore, a total of 521 
DEGs were annotated to "biological process" terms, in which 
71.4% and 63.6% of the assignments belonged to "metabolic 
process" and "cellular process", respectively, followed by "cellular 
metabolic process" and "primary metabolic process" (Figure 2). 

To characterize the functional consequences associated with 
drought response, a pathway analysis of the DEGs based on the 
KEGG database was performed. The results indicated that genes 
participating in the biosynthesis of secondary metabolites, plant 
hormone signal transduction, plant-pathogen interaction, amino 
sugar and nucleotide sugar metabolism, and ubiquitin-mediated 
proteolysis were differentially expressed in PEG-treated group 
(Table S2). 

Genes Associated with the Major Functional Group were 
Affected by PEG-induced Drought Stress 

The global gene expression results showed that a great number 
of genes were significandy affected by PEG-induced drought stress. 
We classified the DEGs into three groups. 

Group A was composed of some regulatory proteins, such as 
transcription factors, protein kinase, and other signaling molecules 
(Table 2). Twenty genes encoding protein kinase, including CBL- 
interacting protein kinase (CIPK), threonine protein kinase, 
receptor-like protein kinase (RPK), mitogen-activated protein 
kinase (MAPK) cascade, phosphoenolpyruvate carboxylase kinase, 
and other kinases, were remarkably changed under drought stress. 
The results showed that six genes were up-regulated, and 14 genes 
were down-regulated. Meanwhile, a total of 14 drought-activated 
or repressed transcription factor genes were identified from several 
different families including WRKY, NAG, MYC, TCP, and bZIP. 
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Table 1. Categorization and abundance of tags. 





Summary 




C 


P 


Raw tags 


Total 


7266032 


7196802 




Distinct tags 


445588 


488325 


Clean tags 


Total number 


7026443 


6913120 




Distinct tag numbers 


212338 


211052 


All tag mapping to gene 


Total number 


5869765 


5356098 




Total % of clean tags 


83.54% 


77.48% 




Distinct tag numbers 


110536 


105096 




Distinct tag % of clean tags 


52.06% 


49.8% 


Unambiguous tag mapping to gene 


Total number 


4261924 


3722859 




Total % of clean tags 


60.66% 


53.85% 




Distinct tag numbers 


85870 


81776 




Distinct tag % of clean tags 


40.44% 


38.75% 


All tag-mapped genes 


Number 


39354 


38413 




% of ref genes 


41.48% 


40.49% 


Unambigous tag-mapped genes 


Number 


27207 


26381 




% of ref genes 


28.68% 


27.81% 


Unknown tags 


Total number 


1 1 56678 


1557022 




Total % of clean tags 


1 6.46% 


22.52% 




Distinct tag numbers 


101802 


105956 




Distinct tag % of clean tags 


47.94% 


50.2% 



Clean tags are tags that remained after filtering out dirty tags (low quality tags} from raw data. Distinct tags are different kinds of tags. Unambiguous tags are clean tags 
remaining after the removal of tags mapped to reference sequences from multiple genes. C represents the control group; P represents the PEG-treated group. 
doi:l 0.1 371 /journal.pone.0094277.t001 



Group B was composed of 13 functional protein genes, which 
included genes involved in flavonoid biosynthesis (6 DEGs), LEA 
protein (2 DEGs), small heat shock protein (sHSP) (2 DEGs), 
aquaporin (2 DEGs), and proline transporters (1 DEG) (Table 2). 
Contrary to the expression of regulatory protein genes, most of the 
functional protein genes were up-regulated after 12 h of drought 
treatment. 

Group C was composed of 14 proteins that were annotated as 
low-temperature-induced protein (1 DEG), dehydration-induced 
protein (1 DEG), defensing precursor (1 DEG), resistance protein 
(4 DEG), universal stress protein (1 DEG), and proteins involved in 
chitinase biosynthesis (4 DEGs) (Table 2). Similar to the alteration 
of functional protein genes, the expression of these 14 genes also 
exhibited a general upward trend. 

Discussion 

Plants, which are sessile organisms, try to respond and adapt to 
stress by altering their gene expression patterns when they are 
subjected to adverse environments [7]. Therefore, the genes that 
showed different expressions are probably involved in stress 
response. Shinozaki and Yamaguchi-Shinozaki [10] reported that 
the drought-inducible genes identified can be classified into 
regulatory and functional proteins. The regulatory proteins 
include various transcription factors, protein phosphatases, protein 
kinases, and other signaling molecules that participate in further 
regulation of signal transduction and downstream gene expression. 
The functional proteins are composed of molecules, such as 
chaperones, water channel proteins, LEA proteins, mRNA- 
binding proteins, sugar and proline transporters, detoxification 
enzymes, and various proteases. In the present study, the DEGs 
between the leaves of drought-stressed and unstressed R. soongorica 



In our results, five transcription factors were identified as up- 
regulated, whereas nine transcription factors were repressed. The 
expression profiles exhibited that the majority of regulatory 
proteins were repressed after 1 2 h of drought treatment. 




U H 1 , 1 , 1 , 1 , 1 , 1 , 

>1 >5 >10 >20 >50 >100 

Tag copy number 

Figure 1. Distribution of total clean tag and distinct clean tag 
copy numbers from the libraries C and P. C represents the control 
group; P represents the PEG-treated group. 
doi:1 0.1 371 /journal.pone.0094277.g001 
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Figure 2. Gene classification based on gene ontology (GO) for DEGs in C and P libraries of R. soongorica. Only the biological processes 
were used for GO analysis. C represents the control group; P represents the PEG-treated group. 
doi:1 0.1 371/journal.pone.0094277.g002 



were investigated to conduct a global analysis of the transcriptomic 
response, which can facilitate our understanding of drought 
tolerance mechanisms. In our results, 1325 DEGs with 379 up- 
regulated and 946 down-regulated genes were detected. Our 
results are similar to that of Fan et al. [6] , who also reported that 
down-regulated genes were more abundant than up-regulated in 
soybean leaves after drought exposure. Functional annotation 
analysis results indicated that many drought-inducible genes with 
various functions were identified in R. soongorica, including a 
number of regulatory proteins, functional proteins, and other 
inducible proteins in response to drought stress. 

Regulatory proteins are involved in signal transduction path- 
ways and in controlling the expression of functional genes in stress 
responses. In higher plants, many regulatory genes, such as those 
encoding protein kinases and transcription factors, are induced by 
environmental signals or stresses [35]. In this study, 20 DEGs 
encoding protein kinases were detected after drought treatment. 
Various protein kinases including CIPK protein, threonine protein 
kinase, MAPK, calcium-dependent protein kinases, and RPK are 
important regulatory proteins in the regulation of some stress- 
inducible genes [36]. Several experiments have demonstrated that 
protein kinase signaling conferred plant resistance to stress by 



regulating downstream functional genes [37,38]. In addition, 
protein kinase also functions in regulating transcription factors 
[39]. A total of 14 genes encoding transcription factors, including 
WRKY, NAG, MYG, TGP, and basic region leucine zipper 
(bZIP), were also significantly difiFerentially expressed. Transcrip- 
tion factors are important in regulating plant responses to biotic 
and abiotic stresses [40]. Altering the expression of certain 
transcription factors can greatly change the expression of a 
number of stress-related target genes that influence plant stress 
tolerance in transgenic plants [9,10]. Thus, the alteration of 
protein kinases and transcription factors may function in the signal 
transduction process by regtilating the expression of various 
functional genes in response to drought stress. 

Aside from signaling transduction, plants adapt to adverse 
environment through the activation or up-regulation of genes that 
directly function in stress resistance [7]. Flavonoids are a class of 
secondary metabolites with diverse functions in growth, develop- 
ment, reproduction, and are also involved in diverse stress 
responses in plants. Flavonoids accumulate rapidly and play a 
protective role when plants are exposed to abiotic stress [41]. In 
recent years, several genes encoding important enzymes, such as 
phenylalanine ammonia-lyase, F3H, dihydroflavonol 4-reductase 
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Table 2. Selected DEGs between the C and P libraries. 




Funtional group Unigene ID 


Gene annotation 


TPM-C 


TPM-P 


Log^Ratio (C/P) 


Protein kinases 


Unigene 9485_A 


CBL-interacting protein l<inase 


2.99 


12.87 


2.11 


CL2201. Contigl_3 


CBL-interacting protein l<inase 


25.9 


59.45 


1.20 


CL4525. Contigl_A 


threonine protein l<inase 


4.98 


15.62 


1.65 


Unigene 2338_A 


threonine protein l<inase 


2.85 


7.96 


1.48 


Unigene 14518_A 


serine/threonine-protein l<inase 


78.13 


35.73 


-1.13 


Unigene 32052_3 


receptor-lil<e protein l<inase 


2.56 


0.01 


-8 


Unigene 39297_A 


receptor-lil<e protein l<inase 


2.13 


0.01 


-7.73 


Unigene 52174_A 


receptor-like protein kinase 


6.26 


1.3 


-2.27 


CL2777. Contig2_3 


receptor-like protein kinase 


36.43 


14.9 


-1.29 


CL6483. Contigl_3 


receptor-like protein kinase 


16.94 


7.67 


-1.14 


Unigene 43240_A 


receptor-like protein kinase 


76.71 


38.19 


-1.01 


Unigene 21624_3 


mitogen-activated protein kinase kinase kinase 


7.26 


1.59 


-2.19 


CL1583.Contig2_A 


mitogen-activated protein kinase kinase kinase 


7.26 


1.59 


-2.19 


Unigene 9581_3 


mitogen-activated protein kinase 2 


25.62 


10.7 


-1.26 


Unigene 633_3 


AlVlP-activated protein kinase 


24.76 


10.7 


-1.21 


Unigene 21618_3 


calcineurin B-like protein-interacting protein kinase 


25.33 


12.01 


-1.08 


Unigene 10286_A 


phosphoenolpyruvate carboxylase kinase 


0.57 


4.63 


3.02 


Unigene 27156_3 


pyruvate, phosphate dikinase 


35.44 


85.2 


1.27 


CL4464. Contig2_3 


pyruvate kinase 


85.82 


29.36 


-1.55 


Unigene 50230_A 


histidine kinase 


40.28 


17.94 


-1.67 


Transcription factors 


Unigene 30342_3 


WRKY transcription factor 7 


4.7 


14.75 


1.65 


Unigene 10063_3 


WRKY transcription factor 6-1 


1.99 


0.01 


-7.64 


Unigene 27699_3 


IVIYC protein 


14.23 


47.74 


1.75 


CL7944. Contigl_3 


NAC domain-containing protein 


14.23 


4.48 


-1.67 


Unigene 31521_3 


NAC domain-containing protein 


12.52 


4.19 


-1.58 


Unigene 25674_3 


transcription factor TCP20-like 


13.66 


33.7 


1.30 


CL7341. Contigl_3 


transcription factor TCP20-like 


45.26 


16.92 


-1.42 


Unigene 22142_3 


transcription factor TCP7-like 


31.88 


15.04 


-1.08 


Unigene 21994_3 


DNA binding protein 


6.55 


14.32 


1.13 


Unigene 23673_3 


transcription factor bZIP124 


51.24 


20.83 


-1.29 


Unigene 1336_3 


transcription factor bZIP12 


4.7 


11.57 


1.29 


Unigene 28584_3 


transcription factor UNE12 


37.86 


12.01 


-1.66 


CL8097. Contigl_A 


transcription factor UNE12 


8.11 


3.18 


-1.35 


CL5168. Contigl_3 


transcription factor UNE12 


12.95 


5.35 


-1.28 


Flavonoids biosynthesis 


CL93. Contigl_A 


UDP-glucose: glucosyltransferase 


22.06 


50.19 


1.19 


CL1676. Contigl_A 


Flavonol sulfotransferase-like protein 


3.56 


9.4 


1.40 


CL2370. Contig2_3 


flavonoid 3'-hydroxylase 


10.39 


25.6 


1.30 


Unigene 1473_3 


flavonoid 3'-hydroxylase 


84.96 


34.72 


-1.29 


Unigene 19573_A 


dihydroflavonol 4-reductase 


1.14 


17.65 


3.95 


LEA proteins 


CL4558. Contigl_3 


late embryogenesis abundant protein 


5.84 


12.73 


1.12 


Unigene 442_3 


late embryogenesis abundant protein group 9 


7.54 


22.57 


1.58 


Heat shocl< proteins 


Unigene 18824_A 


mitochondrial small heat shock protein 


3.84 


10.99 


1.51 


CL12936. Contig1_A 


heat shock protein 17.5 cytosolic class II 


1.85 


8.25 


2.15 


aquaporin 


Unigene 39258_A 


aquaporin 


1.71 


9.69 


2.50 
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Table 2. Cont. 





Funtional group Unigene ID 


Gene annotation 


TPM-C 


TPM-P 


Log^Ratio (C/P) 


Unigene 19242_A 


aquaporin 


426.96 


176.77 


-1.27 


Proline transporter 


CL1109. Contig2_A 


proline transporter 


7.54 


27.34 


1.86 


Cold shock protein 


Unigene 233_3 


cold shock protein 


11.53 


58.44 


2.34 


low temperature-induced protein 


Unigene 24454_A 


low temperature-induced protein 


0.01 


2.03 


7.67 


Dehydration-induced proteins 


Unigene 922_A 


dehydration-induced proteins 


10.67 


21.99 


1.04 


universal stress protein 


Unigene 25385_3 


universal stress protein 


93.22 


263.56 


1.50 


Resistance protein 


CL8753. Contig4_3 


disease resistance protein 


1.85 


7.23 


1.97 


CL6875. Contig2_3 


disease resistance protein 


0.01 


2.03 


7.67 


CLSnS. Contigl_3 


disease resistance protein 


4.7 


1.01 


-2.22 


Unigene 26781 _3 


natural resistance-associated macrophage protein 


1.28 


16.49 


3.69 


Unigene 28389_3 


nematode resistance protein 


36.86 


77.68 


1.08 


defensin precursor 


Unigene 26068_3 


defensin precursor 


20.92 


87.08 


2.06 


chitinase 


CL630. Contig2_3 


class IV chitinase 


3.7 


67.84 


4.20 


CL5695. Contigl_A 


chitinase 


6.4 


41.8 


2.71 


Unigene 30551_3 


class IV chitinase 


13.52 


6.08 


-1.15 


CL630. Contig3_3 


class IV chitinase 


60.91 


27.63 


-1.14 



C represents the control group; P represents the PEG-treated group. The DEGs was selected through FDR <0.001 and |log2Ratio| >1. 
doi:l 0.1 371 /journal.pone.0094277.t002 



(DFR), and atithocyanin synthase, which are involved in the 
flavonoid biosynthetic pathway, were found to be up-regulated in 
response to different kinds of stresses, such as salinity, UV-B 
irradiation, water deficit, and nitrogen stress [41-44]. Increased 
expression of genes involved in the flavonoid pathway was 
detected in our previous study on R. soongorica [16]. In the present 
study, most DEGs encoding key enzymes involved in flavonoid 
biosynthesis were also up-regulated, especially DEGs that encode 
UFGT, DFR, and F3'H. These findings are consistent with the 
results of other studies conducted on birch, potato, and rice after 
stress exposure [45-47]. The up-regulated genes may result in 
accumulation of flavonoids. Flavonoids may serve as signaling 
molecules in signaling cascades. Flavonoids also affect cellular 
function and modulate gene expression in response to adverse 
conditions [48]. Agati et al. [49] found that flavonoids directly 
scavenged the ROS brought about by stress. Therefore, up- 
regulation of flavonoid pathway genes may enhance the protective 
role oi R. soongorica, thus leading to better drought tolerance. 

Secondary metabolites and stress-responsive proteins, such as 
LEA proteins, accumulate naturally in the seeds of some 
desiccation-tolerant plants, and are also induced in vegetative 
tissues when subjected to water-limited conditions [50,51]. LEA 
proteins, largely attributed to their extreme hydrophUic nature, 
have also been implicated in detoxification, sequestration of ions, 
maintenance of protein or membrane structure, binding water, 
and operation as molecular chaperones during dehydration [52]. 
In this study, two DEGs were identified as up-regulated LEA 



proteins. Our results are in agreement with the findings of Gechev 
et al. [53], who found high expression of LEA genes during 
drought and desiccation in leaves of Haherka rhodopensis. We 
concluded that dehydration-induced expression of LEA transcripts 
is a common response in resurrection and desiccation sensitive 
plants [54] . Improved resistance to drought and salinity was also 
observed in transgenic rice when LEA was overexpressed [55]. 
Expression of LEA proteins in yeast confers improved resistance to 
high salinity and freezing [56]. These studies demonstrated that 
LEA proteins functions in binding water, maintenance of protein 
or inembrane structure, or operation as molecular chaperones to 
enhance tolerance to dehydration, whereas few additional clues to 
the mechanism of protein function are provided. 

Plant HSPs facilitate protein folding or assembly under diverse 
developmental and adverse environmental conditions. In this 
study, two small HSPs (sHSPs) encoding DEGs showed up- 
regulation under drought stress. One sHSP was annotated to 
mitochondrial sHSP (MT-sHSP), and the other was annotated as 
cytosolic sHSP17.5. Although the precise functional mechanism of 
sHSPs is still unclear, in vivo studies have shown that sHSPs 
function as molecular chaperones in stressful conditions and in 
normal development [57]. Overproduction of sHSP increased 
drought tolerance in transgenic rice seedlings [58]. MT-sHSP 
protects the electron transport complex I and prevents damage by 
ROS [59] . Lee et al. [60] found that transgenic plants with over- 
expression of MT-sHSP efficiently maintained membrane stability 
by scavenging ROS, and were conferred with enhanced tolerance 
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to salinity. Ectopic expression of cytosolic sHSP 17.1 in Arabidopsis 
thaliana exhibited higher resistance to drought [61]. In a study 
conducted by Sato et al. [62], sHSP 17.5 showed molecular 
chaperone activity in vitro and represented an example of an HSP 
capable of protecting cells against thermal extremes. Therefore, 
these results indicated that the sHSPs may play roles in scavenging 
ROS or functioning as molecular chaperone to respond to drought 
stress. 

The DEGs encoding aquaporins and proline transporter were 
also identified. Aquaporins are channel proteins of intracellular 

and plasma membranes that function in transporting water, small 
solutes (urea, boric acid, and silicic acid), and gases (ammonia and 
carbon dioxide). A specific role for aquaporins in embohsm 
refilling and recovery of stem axial conductance after drought was 
proposed in grapevine [63]. In our results, two aquaporin genes 
were identified. One aquaporin gene was up-regulated, whereas 
the other gene was down-regulated. An overall tendency of down 
regulation for arjuaporin gene was observed. Ilowc'ver, up 
regulation of certain aquaporin transcripts was also detected in 
rice and Arabidopsis leaves [64,65]. Therefore, transcriptional 
control of arjuaporins in drought-stressed leaves appears to be 
more complex [66]. Interestingly, transport of proline may also be 
affected by water deficit. A gene encoding proline transporter was 
up-regulated under drought stress in the present study, which is in 
agreement with the result of Rentsch et al. [67], who also observed 
the induction of a specific prohne transporter in Arabidopsis under 
water deficit. The DEGs encoding aquaporins and proline 
transporter in R. soongorica may ser\'e as an adaptive strategy by- 
regulating water and proline distribution throughout the whole 
plant. 

Moreover, in addition to two dehydration-induced 19 homo- 
logs, proteins that can be induced by other stresses were also 
detected in the study and annotated to low-temperature-induced 
(LT) protein, cold shock protein (CSP), defensin precursor, 
chitinase, disease resistance protein, natural resistance-associated 
macrophage protein, nematode resistance protein, and universal 
stress protein. LT proteins and CSPs are commonly induced by 
cold stress. Various genes are induced by both drought and cold 
stress, suggesting the existence of crosstalk between drought and 
cold-stress signaling pathways [10]. Defensins, chitinase, disease 
resistance proteins, natural resistance-associated macrophage 
protein, and nematode resistance protein were all found to be 
involved in the defense against phytopathogens [68-72]. There- 
fore, the identification of the genes in this study suggested that a 
crosstalk between drought and pathogen attack. The universal 
stress proteins can confer the ability to respond and adapt to 
environmental changes [73]. The identification of all these genes 
suggests crosstalk between responsive genes or pathways and 
multiple abiotic or even biotic stresses [23]. Moreover, these 
identified genes may provide gene resources for understanding the 
drought tolerance mechanism of iJ. soongorica. 

Among the identified DEGs, most regulatory proteins were 
down-regulated, whereas the majority of functional proteins were 
up-regulated. The different expression patterns shown by the two 
groups are interesting and worthy to be explored. Therefore, we 
tried to discuss the probable interaction of regulatory proteins and 
functional proteins. Most gene regulators exhibit either antecedent 
or simultaneous change in the expression level when compared 
with their targets [74]. One reason is that it takes time for the 
regulator gene to express its protein product and affect (directly or 
indirectly) the transcription of its target gene [75] . Other reasons 
may have been caused by the different mRNA half-lives of the 



regulator and target genes, and the different profiles between 
regulator genes and functional genes may also indicate a feedback 
loop, where a target gene can regulate its regulator [76]. Given 
that signaling pathways are complex dynamic events that occur 
over time, single time point expression profiles in our study are 
insufficient to elucidate temporal events. The different expression 
profiles of regulatory and functional proteins observed in this study 
are interesting. Expression profiling experiment with a series of 
time points should be performed to elucidate this problem in 
further studies. 

Conclusions 

Based on the putative function analysis of annotated DEGs, 
DEGs are classified into three groups. The first group, which 
includes protein kinases (CIPK and threonine protein kinase) and 
transcription factors (WRKY, MYC and TCP), regulates signal 
transduction. The second group includes genes involved in the 
biosynthesis of flavonoids in scavenging ROS, sHSP, LEA 
proteins, aquaporins, and proline transporter that represents 
functional genes that directly enhance drought-stress tolerance. 
The sHSP maintains protein or membrane structure and act as 
molecular chaperones that protect cells from injury under drought 
stress. Aquaporins and proline transporter regulates water or 
proline distribution throughout the whole plant. The third group 
contains proteins induced by other stress factors, which suggests 
that a crosstalk between different stresses. In conclusion, the 
tolerant ability of R. soongorica may hv gain(^d through effective 
signal transduction and enhanced protection of functional proteins 
to reestablish cellular homeostasis. 

Supporting information 

Figure SI Sequencing saturation analysis of the two libraries. C 
represents the control group; P represents the PEG-treated group. 

(TIF) 

Figure S2 Mapping of distinct clean tags in the two DGE 
libraries. C represents the control group; P represents the PEG- 
treated group. 
(TIF) 

Table SI Gene classification based on gene ontology (GO) for 
DEGs in C and P libraries otR. soongorica. C represents the control 
group; P represents the PEG-treated group. 

(XLSX) 

Table S2 KEGG analysis for DEGs in C and P libraries of R. 
soongorica. C represents the control group; P represents the PEG- 
treated group. 
(XLSX) 
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